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We present in this paper an analytical model for a cold bosonic gas on an optical lattice (with 
densities of the order of 1 particle per site) targeting the critical regime of the Bose - Einstein 
Condensate superfluid - Mott insulator transition. We focus on the computation of the one - body 
density matrix and its Fourier transform, the momentum distribution which is directly obtainable 
from 'time of flight" measurements. The expected number of particles with zero momentum may be 
identifled with the condensate population, if it is close to the total number of particles. Our main 
result is an analytic expression for this observable, interpolating between the known results valid for 
the two regimes separately: the standard Bogoliubov approximation valid in the superfluid regime 
and the strong coupling perturbation theory valid in the Mott regime. 

PACS numbers: 



I. INTRODUCTION 

Since their experimental realization in 1995 1], Bose - Einstein condensates (BEC) have become one of 
the most exciting fields in physics. Because the high degree of control and the good understanding of the 
microscopic physics involved, they provide an excellent opportunity to investigate various issues in atomic 
and molecular physics, quantum optics, solid state physics and even high energy physics and cosmology 0]- 

The interest in these systems is also boosted by its possible use in the implementation of quantum in- 
formation processing (QIP)|3|- Cold neutral atoms in optical lattices are a naturally scalable system, and 
because of the weak coupling to the environment long decoherence times are expected. There are detailed 
proposals on how to build quantum gates M ^'^^ qubit buses 0] to exchange information between different 
locations. All these properties make these systems a promising candidate for QIP. 

In most proposals, the physical qubit is a single atom which may be in one of two preferred hyperfine 
states. This implies a strict control of the number of atoms per site, which in principle may be achieved by 
driving the system deep into the Mott insulator (MI) regime 01 ■ However, the gas is usually first condensed 
in a trap, and then the lattice is imprinted on it. This implies driving the system through the superfluid (SF) 
- insulator transition. As with other phase transitions, we expect the particle distribution will be determined 
by events at or just below the critical point; once the hopping parameter is low enough, this distribution will 
be simply frozen in Q . 

To amplify this important point, we observe that it is expected both Landau and Beliaev damping will 
be strongly suppressed in the Mott regime Q; this means that the equilibration times will grow sharply as 
we cross from the superfluid to the insulator phases. The pattern of correlations among different sites and 
particle number fluctuations will get frozen once the relaxation time is long compared with the characteristic 
time in which the parameters of the model are being changed. Unless this change is made very slowly, this 
will happen soon after entering the Mott regime. In this "diabatic" transition, the likelihood of a vacancy 
or of a multiply occupied site will correspond to those of a lattice near the critical point, rather than to the 
parameters of the operating regime. 

The goal of this paper is to formulate an analytical model for a cold bosonic gas on an optical lattice 
(with densities of the order of 1 particle per site) targeting the critical regime of the BEC superfluid - 
Mott insulator transition 0, We focus on the computation of the one - body density matrix ^2 
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and its Fourier transform, the momentum distribution which is directly obtainable from 'time of flight" 
measurements 

m 111 111 (see Hi). The expected numbe r of particles with zero momentum may be 
identified with the condensate population, if it is close to the total number of particles. Our main result is 
an analytic expression for this observable, interpolating between the known results valid for the two regimes 
separately: the standard Bogoliubov approximation valid in the superfluid regime and the strong coupling 
perturbation theory valid in the Mott regime [17i. .18t . .19,. 201 ■ Comparison of our analytic results with exact 
numerical solutions for N particles in a one-dimensional lattice of = 9 sites shows that unlike the standard 
Bogoliubov and strong coupling perturbation our analytic solution sustains an uniform accuracy throughout. 



A. The model 



We consider a system of TV particles distributed over lattice sites, with an integer mean occupation 
number n = N/Ng. In terms of the creation and destruction operators (t) and (t) , the dynamics is 
described by the Bose - Hubbard Hamiltonian (BHM) 



where the first term describes hopping between sites, and the second term the in-site repulsion between 
particles. The matrix Jij is equal to J if the sites i and j are nearest neighbors, and zero otherwise. When 
the repulsion term dominates, the ground state of the system has definite occupation numbers for each site, 
and weak correlations among different sites. The system is in the so-called Mott insulator (MI) phase. When 
the hopping term dominates, atoms condense into a single quantum state extended over the whole lattice; 
the system is in the superfluid (SF) phase. 

In this paper we shall focus on the calculation of the one - body density matrix 



and its Fourier transform 



cTiik ^ (4 (t) ak (t)) , (2) 



Ik 

Nq is the expected total number of particles with momentum q (in units of h/NsO., where a is the lattice 
spacing). Nq may be identified with the condensate population. 

In the deep Mott regime (J = 0), auk = n5ik and Nq = n is the same for all modes. In the opposite limit 
{U = 0) auk = n for every pair of sites and Nq = N6qo. 

Our goal is to obtain analytic expressions for this observable in the intermediate regime U/nJ ^ 1, with 
n ~ 1 as well. 



B. Some approaches to the one-body density matrix 



To motivate our perspective below, let us begin with a brief discussion of some of the most common 
approaches to this problem in the literature and place our work in this context. We feel that other than the 
few full-fledged numerical calculations ([l^llS^), none of the analytic approaches fully cover the transition 
regime described above. Moreover, even if a numerical calculation is feasible it is useful to have a reliable 
analytic approach to match against. 
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To begin with, since our interest is ui, approaches based on the GutzwiUer ansatz or mean field theory 
p^l would not be sufficient. These methods are very powerful to investigate the phase diagram, but because 
they treat different sites as independent, they severely distort the one-body density matrix. 

These approaches may be improved on, of course. The GutzwiUer ansatz may be taken as just the first 
step in a consistent perturbative expansion |25| , and the mean field decoupling ansatz may be applied to full 
cells rather than individual sites (26| . However, the required order in perturbation theory (or the size of the 
fundamental cell) to get a reliable result scales with the size of the lattice, and soon the difficulty becomes 
comparable to a full numerical solution. 

Starting from the superfluid regime, the simplest way to get cri is the Bogoliubov approach ^27j . Since we 
shall consider the case in which the gas is at fixed total particle number, rather than fixed chemical potential, 
we must consider instead the particle number conserving (PNC) formalism '28l|. However, for the purpose 
of this preliminary discussion we may make abstraction of the difference. 

A simple minded mean field approach, in which we simply replace aj by its "expectation value" Zj, is 



dobal phase invariance aj 
|. In other words, simple- 



bound to fail. Since the BH Hamiltonian has the ; 
theorem the mean field theory must be gapless |2£ 
only describe the superfiuid phase. 

Since the one - body density matrix is the time coincidence limit of the two - 

one could think of finding equations of motion for these functions directly, without including a mean field 
[30l| . but this approach also fails. In a nutshell, the difficulty is as follows. The Heisenberg equation of motion 
for a| (t) is 



e Qj, in view of Goldstone 
minded mean field theory can 

point function (t) (f) 



H,a] (t) 



(4) 



whereby 



dt 



(a] it) a, it')) = (^')) - u{{^f^o) it) it')) , 



(5) 



and we face a closure problem, namely, how to express the four point function in the last term in terms of 
two point functions. A typical resolution is a Hartree - like scheme, where we approximate 



[[afa,'^ (t) ak [t')) ^ 2 (a] {t) a, (<)) (a] (t) a,, (t')) =^ 2n (a] (t) (i')) ■ (6) 
However, in the weak hopping limit we expect the system will be close to the MI ground state 



\MI)=l[\n).^, 

i 

(that is, each site is in a state of weh defined occupation number) where we can compute 



(7) 



a] (t) ak (t')) ^ "n^i. 



(8) 



[(afa,^ (t) ak {t')) ^n{n^l) % « (n - 1) (a] (t) ak {t')) . 
We see the Hartree approximation is off by a factor of two, even if n ^ 1 [sJl ■ 



(9) 
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A possible way around this problem is to obtain a formal equation of motion for an object (say a two 
point function) for finite U and J, and then approximate the coefficients in the formal equation (for example, 
a self-energy) by their exact value at J = or for very large J, as needed 0, [sj- However, the actual 
expressions derived in this way are not reliable at the transition region, which is where our main interest 
rests. 

In the opposite Mott insulator regime, the most straightforward approach is Rayleigh - Schrodinger per- 
turbation theory in the parameter J [l7l flsL Il9ll20| . However, the complexity of the calculation increases 
steeply with each increasing order, and so its accuracy for finite values of J is hard to assess. Comparison 
against exact solutions for n = 1 and Ns = 5, 7 and 9 shows that first order perturbation theory breaks 
down before the transition (see below). This is consistent with the expectation that perturbation theory 
breaks down when Jn > U. 

Dilute gases with very strong repulsion may be treated as a free Fermi gas |3^ . This approach has been 
recently successfully extended to densities n > 1 [33 |. 

Returning to our the above failed Hartree attempt, it is clear that the closure problem arises in the U 
term because it is the nonlinear term, while the J term is linear. One obvious alternative is to reformulate 
the theory in such a way that this situation is reversed. This is accomplished in the so-called slave boson 
/slave fermion method (35j . 

The slave boson method requires the introduction of a large number of auxiliary fields and new constraints 
on the theory. In this paper we shall explore a similar strategy (that is, making the repulsion term linear, 
the hopping term nonlinear) while keeping closer to the original fields in the Hamiltonian. 

One possible way to implement this is to observe that the interaction term is actually quadratic on the 
site occupation number Ui = aj^a;, since a^af = rii {rii — 1) . This suggests to consider as fundamental a 
"phase" variable cpi canonically conjugated to the occupation numbers rii HE US 

[nj,ipi] = -iSij, (10) 
(here and after we assume ?j. = 1). The original creation and destruction operators are 

aj = [exp-i(pt] (11) 



al = [expzc/?j] . (12) 

The implementation of this idea hits some well known difficulties |38) . If the operator (pi exists and is 
hermitean, then the operators exp—ir(pi are unitary and shift the state |n) into \n — r). But such operators 
annihilate the vacuum state |0) , so they cannot be unitary. We shall return to these difficulties below. 

In terms of the density and phase variables, the classical Hamiltonian becomes 



if = ^ < - ^ 2Jijy/n~n~ cos [ipi - ipj] + [n^ - 1) Z • (13) 

If we further approximate ^Jn^nj = n in the hopping term, then we obtain the quantum phase model 
|39l l40l |. This model displays a phase transition, and it has been used to investigate noncquilibrium aspects 
of the Mott transition ||8|. 

On closer examination, the approximation involved is valid when Un > J |4l|. Therefore, for n ^ 1 it 
fails at the transition region. In conclusion, while the quantum phase model is the best option on the shelf, 
it must be generalized to lower densities to be truly reliable in the relevant regime [i^. 

One possibility is to allow for particle fluctuations, but only as far as any given site is never more than 
one particle above or below the average. Then it is possible to map t he p roblem onto the XY model or else 
use a path integral representation in terms of spin 1 coherent states l4!^ . These model also display a phase 
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transition, and a Gross - Pitaievsky description has been recently developed. However, we are not aware of 
attempts to carry the perturbative evaluation of these models to higher orders. Below we shall explore an 
alternative strategy with the same overall goals. 

Finally we observe that the so-called truncated Wi gne r approximation and other phase space methods 
have been successfully applied in the n >> 1 limit 44, 45]. 

From this description we see the lack of suitable treatment in the literature of the one body density 
matrix at the transition region for low densities. Not only there is no single approach which is fully reliable 
throughout, but moreover those which are successful on one asymptotic regime are based on a quite different 
physical model than the ones which succeed on the other (compare, e.g., Bogoliubov methods against the 
Tonks - Girardeau gas approach or strong coupling perturbation theory) . A model which is able to describe 
the transition region within a single physical model and keeping an uniform accuracy would be a definite step 
forward. This is our aim here. To be fully understood, however, we must identify some desirable features 
any new approach to this problem must possess to be truly useful. 

C. Our approach in the context of ongoing research 

As we mentioned above, our interest in this problem of the loading of BEG atoms on to an optical lattice 
is motivated by the feasibility of using this process to initialize a quantum computer. This longer range goal 
sets certain constraints on the model which we choose to perform our analysis. 

The first consideration is that, although in this paper we shall only discuss the equilibrium case, in last 
analysis one needs a full nonequilibrium formulation of the problem. With this goal in mind, we adopt the 
Schwinger - Keldysh or closed time - path (CTP) 30, 46] formalism from scratch. As a side benefit, we 
shall see below that this choice is also helpful in overcoming the formal difficulties of the density - phase 
representation. 

A related requirement is that there should be a well defined way to carry the perturbative evaluation of the 
model to any order, but because this will be unavoidably complex, already the first order in the expansion 
must give sensible results. In particular, it is desirable to have the model in path integral language, as it is 
most adapt to further implementation of perturbation theory. 

Actually, the simplest quantum phase model formulation fails this test; with some oversimplification, the 
problem is that \/n + 5n ^ ^/n + Sn/2y/n is a bad approximation ii 5n > n 47]. We shall seek a new set 
of variables in which the perturbative evaluation of ai is more accurate than in the original ones. We shall 
show this by comparing the first order approximation to our model with the exact solutions in the case of 
small systems, and to the PNG and strong coupling perturbation theories for larger densities. 

It is seen in actual experiments that collisions with noncondensed particles and loss are not significant 
except on the longest time scales (above Is ^3)- Therefore we shall consider the case of an isolated gas, 
i. e., the total number of particles will be constant as opposed to the case of a gas interacting with a 
particle reservoir, whereby the chemical potential remains constant. However, instead of the PNG approach, 
we shall develop a formalism which is more suitable to the path integral formulation of the model. We shall 
regard the given value of the total particle number N = nNg as a constraint on allowed states of the system, 
rather than just a dynamical condition. The resulting theory will amount to an independent quantization 
of the system; our model and the PNG one will agree only with respect to the time evolution of observables 
which commute with N. Of course, a| (t) Ui (t) is one of these observables; not so the creation or destruction 
operators separately. A detailed comparison of the path integral and PNG approaches is given in Ref. j43 . 

Let us observe that this procedure is less unusual than it may seem. For example, in studies of the ground 
state of the system, it is common to adopt trial wave functions which preclude site occupations farther from 
the mean than a few units (a similar policy is sometimes adopted for the numerical diagonalization of the 
Hamiltonian) . In practice, this means that the Hilbert space of allowable states is constrained; the reduction 
is actually more drastic than the one postulated here. 

A similar procedure has been implemented in the field of high temperature superconductors near the Mott 
limit, to enforce such constraints as excluding double occupancy [SOj ]. 

From the technical point of view, the advantage of taking the given value of iV as a hard constraint is that 
in the constrained system, the global phase invariance of the Bose - Hubbard model (BHM) becomes local 
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in time. Technically, the model becomes a gauge theory, with the constraint N as gauge generator |5l|. This 
allows us to take advantage of the powerful methods of gauge theory quantization (of which we shall only 
have a glimpse in this first take of the problem) [s^ . 

When seen under the light of our stated long term goal, a number of shortcomings of our present work 
clearly stand up, and it is only fair that we mention some of them. First, to be sure we have an accurate 
description of the transition we should also compute other observables, such as the particle number fluctua- 
tions ^3 and the dynamic structure factor |20ll5^ . We have only considered a ho mog eneous lattice, while 
a lattice superimposed to a harmonic trap would be more relevant to applications 55] (the presence of the 
trap has a drastic effect on the transition H^). We have not considered lattice fluctuations 57] nor finite 
temperature effects [H^. We have considered a condensate of atoms without internal degrees of freedom, 
while of course the internal structure is essential for QIP |5^ . 

It is also interesting to observe that some of the quantities we compute in this paper have been measured, 
both in one and three dimensional systems |6^, '6l']. We will comment briefly on these results in Section VII; 
a detailed discussion will be given elsewhere 62]. 

In spite of these unachieved goals, the formulation of a fully analytic theory of the one body density matrix 
is a necessary first step towards constructing a realistic theory of the loading process, which we now proceed 
to tackle. 



D. This paper 



The rest of the paper is organized as follows. Over the next four Sections, we develop a formal presentation 
of the model. In Section II, we develop the CTP path integral representation for expectation values of BEG 
observables. In Section III, this representation is translated to the density - phase representation. In 
Section IV, we shift to a new set of variables, more suitable for the further perturbative evaluation of these 
expectation values. In Section V we explore the simplest approximation, where the theory is linearized in 
the inhomogeneous modes. 

In Section VI we apply this machinery to the computation of the one - body density matrix and the 
momentum distribution function. In the final Section VII we show the results of a comparison of om model 
against both exact numerical results and other approximated approaches, and conclude with some final 
remarks. 

In Appendix A we present a brief derivation of the other approximate approaches discussed in the Results 
Section, namely, first order perturbation theory in J/N and the PNC approach to first order in iV~^. This re- 
sults are not new, and are included only to prevent any misunderstandings due to different notations between 
this and the original papers. Appendix B discusses the validity of Eq. (|142|l below as an approximation to 

Eq. (innj). 



II. THE CTP PATH INTEGRAL APPROACH TO BECS 

In this Section we'll put together the basic formulae for the coherent-state path integral method [g^] to 
compute expectation values of observables within the causal CTP approach |4q ]. 

Before we get down to formulae, let us try and convey the idea of the approach in simple terms. Let 

us begin with the problem of computing the vacuum expectation value of some observable O. One 

possibility is to add a new term — J O to the Hamiltonian. Let us call H the original Hamiltonian, Hj the 
new Hamiltonian Hj = H — J O. Then, if we can find the ground state energy Ej of the new hamiltonian, 

first order perturbation theory impliest that = —dEj/dJ at J = 0. We have translated the problem 

of computing into the problem of computing Ej. 

As it turns out, a surprisingly efficient way of computing Ej is by computing the matrix elements of the 
euclidean evolution operator |64| 
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= (0 



-,-TH, 



0) 



(14) 



where |0) is the vacuum state (we may assume that the external source J is swithched olF adiabatically at 
infinity, so the vacuum is unambiguously defined). It turns out that when the euclidean lapse T — s- oo, 

We [J] -TEj, so again (^O^ = T^^dWe [J] /dJ at J = 0. 

In a time-dependent situation, however, an euclidean formulation is not readily available. One may attempt 
to make do with the analytical continuation of Eq. H14|) back to physical time {T inside the bracket is the 
time-ordering operator) 



't-^l = (OOUT 



T 



^-i J dt Hj 



(15) 



but this is untenable. The "expectation values" derived from Wm are generally complex, even if O is 
Hermitian, and they do not evolve causally in time [65l |. The problem comes from the fact that they are 
no longer expectation values, but rather matrix elements between the asymptotic vacua |0/iV) and \0OUT), 
which are not necessarily equivalent in this time-dependent situation. 

The solution found by Schwinger (46l | was not to include one external source but two, and J^, and to 
define a new generating functional 



^iWcTp[J] _ 



= (OIN 



T 



i J dt Hj2 



T 



i J dt H ji 



(16) 



where T is the anti time-ordering operator. We may also think of and as a single source defined on 
a "closed time-path" which reaches from t = 0, say, to the far future (wherein the source takes the values 
{t)) and then bounces back to i = (the source swithching to (t)), therein the name of the method. 
It is readily shown that differentiation of Wctp yields true expectation values, which are of course real and 
evolve causally. 

Since both quantum states in Eq. (I16|l are defined at the same reference time i = 0, Wctp is readily 
generalized to non vacuum situations 66] . Let pi be the density matrix describing the state at t = ti. Then 
expectation values with respect to pi may be obtained from the CTP generating functional 



iW 



ti{U^^ {tf,U)Ui{tf,U)p{U)} 



where 



U^{tf,U 



T 



— i r'^' dt H'{t) 



(17) 



(18) 



Our problem is to build a generating functional which will allow us to compute the one - body density 
matrix. Our starting point will be Eq. (|17|l . For reasons of efhciency, we shall seek a path integral 
representation of the trace. 

In this and the next two Sections, we shall construct this representation. In this Section we shall use 
the well known coherent state representation |63l |. putting the emphasis on the implementation of CTP 
boundary conditions. Then we shall proceed to rewrite the CTP generating functional in terms of more 
suitable variables, to optimize the accuracy of its perturbative expansion. 



A. The coherent state representation 



We shall begin by recalling the usual coherent state path integral representation of transition amplitudes 
[63l| . The CTP boundary conditions shall be introduced in next Subsection. 
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For simplicity, let us consider a single one-particle state. There is a basis made of occupation number 
eigenstates \n) 

N\n)=n\n) (19) 
In particular, there is the vacuum state |0) . These states are orthonormal and complete 

(m \n) = 6mn (20) 

J2\n){n\ = l (21) 
The destruction and creation operators relate states of different occupation numbers 

a\n) = s/n\n-l); \n) = \/n + 1 |n + 1) (22) 

Therefore 

a^a = N; [a,at] = 1 (23) 
A coherent state \a) is an eigenstate of the destruction operator 

a\a) =a \a) (24) 

Adopting the normalization (0 |a) = 1 we find 

{n \a} = — (25) 

Let \b) be a second coherent state; then 



{a\b) =exp{a*b} (26) 

The vacuum is the coherent state with a = 0. 

While not orthogonal, the coherent states are complete, in the sense that 

/ ^exp{-a*a}l«>H = l (27) 
We may use the completeness relationship to write down the trace of an operator A 

tvA = y (n| A |n) = / ^^exp {-a*a} {a\ A \a) (28) 
^-^ J zm 

Now consider the transition amplitude between the state |aj) at time ti = Q and the state |a/) at time tf, 
where | a/) is the eigenstate of the Heisenberg operator a{tf) with proper value a. Since a{tf) = e^^^f ae~^^*f , 
we have {h = 1) 
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\-af) = e^"'f \a) (29) 

and 

{af\a,) = {a\e-'"'f\a,) (30) 

Let N be some large number and e = tf/N. Write at — ao, a — un- Then, inserting — 1 identity 
operators, we have 

(a/ = 7 I n ^^exp{-a>„}(a„+i|e-'^^|a„)| (ai|e-'^^^|ao) (31) 
which may be written as (assuming the Hamiltonian H = H (a^ , a) is in normal form) 

(5/ |a,) = j [Da]^_^ exp {i^jv [a*, a]} e"«"" (32) 

where 

N-i , * , 

[Da],_, - n ^ (33) 

Tl=l 



Sn [a*, a] = ^ {iflJi {an - fln-i) - ei? (a,* , ftn-i)} 

n=l 

Going to the contimmm limit, where o„ — a„_i ~ edajdt^ we get 

<Ja(t/)^ |o(t,)) = I {Da\ exp{iS'[a*,a]}e''*"(*^) (34) 

(35) 



S'[a*,a] = y" (it |m*|^-iJ(a*,a)| 



The integration is over paths where the initial value of a and the final value of a* are fixed, and given by 
a (ti) and a* {tf) , respectively. 

B. The CTP boundary conditions 

We now have all the necessary elements to evaluate the CTP generating functional Eq. p7|) . The idea 
is that the initial density matrix p is propagated forwards in time with some Hamiltonian and then 
backwards with a Hamiltonian . Insert three identity operators in Eq. p7l) to obtain 

iw f da%daN dah*da}) daQ*dan r / « i* i 9* 9\i 

e- = J ^- 2_Jl^_Jlexp{-(aW + 4*«J + a^*«^)} 

{qnI U2 {tf,U) \aiy {unI Ui (tf.ti) \al) {al\p{ti) \al) (36) 
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Now use the corresponding path integral representations 




(37) 



The configuration on the forward branch has (0) = aj and o^* (tf) = a^. On the backward branch, we 
have a^* (0) = flg* and a? {tf) — a^. Once W is known, causal expectation values may be computed by 
differentiation. 

Eq. (|37|) is the main result of this Section. In order to make use of it, however, we must rewrite it in a 
more suitable set of variables. This translation is the subject of the next two Sections. 



In this Section we present the basic elements of the path integral formulation of a system of bosonic atoms 
in an optical latice in terms of number and phase variables, while enforcing a fixed total particle number. 
This will set the stage for a further canonical transformation to a more convenient set of degrees of freedom, 
to be carried out in the next Section. 



Our starting point is the Madelung representation for the creation and destruction operators, Eqs. (|11|) 
and H12|) . The phase observables ipi have eigenstates \ipi) , which are a complete basis if (pi runs over a full 
circle. To account for the periodic nature of these variables, we define the inner product 



with k running over all integers. A transition element is decomposed into transitions between phase eigen- 
states, mediated by these identity operators. However, as shown by Kleinert 67], all but one of these sums 
may be avoided if we allow tpi to run over all real numbers, and not just a circle. Since the line is the covering 
space of the circle, we shall call this extended theory the covering theory. 

In the covering theory, the discrete observable Ui is replaced by a continuous observable p^, whose eigen- 
states \pi) [such that {p[ \pj) = 6ij6 {p' — p)]) generate the Hilbert space. Then the physical subspace is the 
one generated by the \pi) where pi happens to be a nonnegative integer. 

In the expanded Hilbert space, we have the p representation of a state 



III. DENSITY AND PHASE VARIABLES IN THE CTP FORMULATION 



A. Madelung representation for the creation and destruction operators 




(38) 



k 



{Pi IV') = V'(P») , 



(39) 



In this representation, the operator 



d 



(40) 



opi 



meaning that 



(41) 
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Therefore 

{pi |exp(-zrv5i)| -0) ='4>{pi+r). (42) 
So if IV') — \n) , il) (pi) — d {pi — n) and exp(— irc/j^) \n) — \n~r),d& expected. 

d 

{Pi Wi\ ft) = i-g^ {Pt I fi) = </J» (Pi I fi) , (43) 



(p. |^.)=e-^'^"'-. (44) 

The reason this scheme works better in the CTP formulation than in other approaches is that it affords 
us the freedom to choose an arbitrary initial condition. As long as the initial condition is chosen within the 
physical subspace, the dynamics of the system itself warranties that there will be no unphysical results. In 
particular, we turn the covering theory into the physical theory by inserting the one missing identity of the 
form Eq. into the path integral in a suitable way (see below). 

B. CTP path integrals in number and phase variables 

Before considering the BHM, let us discuss how to build CTP path integrals for the BHM in terms of 
number and phase variables. The key is to clarify the boundary conditions the histories within the path 
integral must satisfy. 

Our starting point is the path integral representation Eq. I|37f) for the CTP generating functional. The 
action is given in Eq. (|35|l . whereby ia* is formally the momentum conjugate to a. To transform the action 
to density and phase variables, define a generating functional 

Q = ^a2p2^v (45) 

so 



Then 



and 



dQ dQ 



ia*da - Hdt = pdip - Hdt - dQ (47) 



S= I dt \p^~Hip,^)}-[Q (tf) - Q iU)] (48) 



Jdt[, 



H = -Y^ Jijy^p7pJ[expi [ipj - ip,]] + ^p, {pi - 1) (49) 
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where Q = {-i/1)N. 

The other factors in the measure may also be expressed in terms of N. In this representation, we consider 
paths which begin at values (0) and (0) and end at a common phase (T) = ip^ (T) = ipf- Observe 
that both the initial and final values of the densities are undetermined. 

Of course, the physical phase variable Lp must be identified with periodicity 27r. We do this by inserting 
the identity operator Eq. (|38|l at some point within the path integral. So we have two kinds of expectation 
values: the expectation values of the covering theory (without the insertion) and the physical expectation 
values (with the insertion), which we shall call (()) and () , respectively. The relationship between these 
constructs will be further clarified below. 



C. Enforcing a fixed particle number 

The quantum theory of the BEG may be regarded as the quantization of the nonrelativistic classical field 
theory defined by the action functional Eq. (|48|l where the canonical variables are pi {t) and their conjugate 
momentum ifi (t) . We are interested in the case in which particle number takes on a definite value N . We 
may reinforce this point by adding a constraint on the theory. This is achieved by introducing a Lagrange 
multiplier /i {t) , and rewriting the action as 



Sftxed = S + J dt p{t)'Y^{p^ - n) (50) 

i 

The original action Eq. (|48|l is invariant under a global transformation tfn [i) —> ipi (t) + constant but the 
new action Eq. (|50|l is invariant under the local (in time) transformations 



>fi^{t)^^,{t)+e{t), (51) 

When is infinitesimal, these are just canonical transformations generated by the constraint. Therefore it 
must be quantized using the methods developed for gauge theories, such as the Fadeev-Popov method. 

This comes about because now the path integral is redundant, since we may transform the fields as in Eq. 
(|51|l . We may fix the redundancy by factoring out the gauge group. Choose some function fg = f [pg, ipig] , 
such that dfg/d0 ^ 0. Then 



jw 



e / Dpt (t) D^t (t) Dp'' 



\Det 



'Sfg 



(52) 



where 



J D9 (53) 
is the volume of the gauge group we wish to factor out, 



S?ot = S^ + J dtY^pit) [ft -n] + -Jdt /2 (54) 

i 

and s is the "gauge fixing parameter" , which may be chosen freely. The determinant may be expressed as a 
path integral over Grassmann "ghost" fields |5^ . For simplicity, we shall adopt a gauge fixing condition / 
which transforms linearly. 
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so that its determinant is a constant and may be ignored. 

In the path integral, /i^ (0) and (0) are integrated over. In principle, there are no restrictions at T, but 
we may assume /i^ (T) — fi^ (T) with no loss of generality. Physically, /i has the meaning of a fluctuating 
chemical potential. 

We note that the freedom to choose the gauge fixing condition / and the gauge fixing parameter s is the 
key to the power of the method. In this paper, we shall restrict ourselves to the simple choice above for 
/ and to the "Landau" gauge s — > 0. Other choices may be used to meet the demands of more advanced 
applications or to optimize the convergence of perturbation theory. 

A different strategy to introduce freely chosen functions in the formalism and then exploit the freedom 
therefrom is the so-called stochastic gauge method 68] . 

Generating functionals as defined in Eqs. H37|) and (|52|) . after adding an external source coupled to an 
operator O, may be used to generate correlation functions involving this operator. If O commutes with 
N both representations are equivalent. Otherwise, they will yield different results. Fortunately, when we 
compute the one body density matrix we are within the domain of equivalence of both formalisms. 



D. Vanishing of the order parameter 

As a check on the formalism being developped, let us verify that the order parameter (a^ (t)) vanishes 
identically. This must hold in any system with a finite number of particles. 
To compute the order parameter, let us add a new source coupled to ipl 



dt' EjK^Ov'KO 



(56) 



Now observe that 



where 



4 (0) = (V7^(0)j, 



iVP^ it))-,, = / ^P" (t) D^t it) (t) exp <^ I 



■^tot ^tot 



dt 



3k,{t')v]{t') 



and 



(57) 



(58) 



jki (t) = -6 {t- ti) Sik 

Let us now show that this vanishes. Make the change of variables within the path integral 



(59) 



a 

^nt) ^ it) + SvUt) : 5]<5^^(i) = (60) 

i 

The homogeneous phases if"" (t) appear linearly in the action. When we integrate them out, they enforce the 
constraints 
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E 



dt 



(61) 



But these constraints are impossible to meet, since they contradict the further constraints from the integration 
over fi in the s — > gauge. 



In this Section, we perform a canonical transformation from the phase and density variables introduced 
above, to a new set of degrees of freedom which are more adept for the perturbative evaluation of the one - 

body density matrix (aj {ti) Uk (^i)) • 



As in the evaluation of the order parameter in the previous Section, the basic idea is to write the creation 
and destruction operators in their polar representation, and then consider the exponential terms as the result 
of external sources coupled to the phases. 

In the general case, a closed evaluation of the one-body density matrix is not possible. We may try assuming 
that all quantities may be decomposed into an homogeneous component and a small inhomogeneous part 



etc. However, one is concerned about the ^/pi factors in the expectation value. Concretely, while the action 
becomes quadratic in the J limit, and so a linearized approximation would seem reasonable at small 
enough hopping, a term like ^/pl does not become Gaussian in any controlled way. 

Our proposal is to introduce a new set of canonically conjugated variables, so that no square roots appear 
in the definition of the one - body density matrix or the Hamiltonian. This will make the perturbative 
expansion starting from a quadratic approximation to the Hamiltonian more straightforward. 

However, as stated in the Introduction, it is not enough to have a well defined perturbative expansion, but 
already the first terms must give sensible results. In our case, the first term in the expansion corresponds 
to keeping only the quadratic terms in the Hamiltonian. This simplified model will be investigated in next 
Section. 



To avoid the nonanalytic square roots in the creation and destruction operators, we proceed as follows. In 
the first branch, we define a new (complex) variable from 



IV. A NEW SET OF DEGREES OF FREEDOM 




pi = pi + 6pl 



(62) 



A. 



A new set of variables 




(63) 




(64) 



This is actually a canonical transformation, since 




,1 ' 




(65) 



The new conjugated momentum is again p;^ 



<l. It follows that on the first branch 
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On the second branch we write instead 



.2t 



Again there is a generating function 



Q^ = ^^[exp2z(e-^?)]=^E^? 

i i 

SO p| is the momentum conjugated to The action 



In the second branch, therefore 



ExpHcitly, the Hamiltonians read 



al = [exp -i£,f*] pf 



H {p\e) = ^ E J^^p^ [-p^ [^1 - m + E U (pI - 1) 



(66) 



(67) 



(68) 



(69) 



(70) 



(71) 



H (p^ e*) = - E [^-p^ [^f - ^'1] p'+H ^Pi (Pi - 1) 



(72) 



plus the gauge terms, in both cases. Observe that in the new variables, the action is explicitly analytical. 



1. Canonical matters 



If we regard the Oj as q-numbers with equal time commutators 



we conclude 



[pi,Pj] = 



(73) 



(74) 



[pi,aj] = -aiSi. 



(75) 
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So 



Now observe that 



but also 



Therefore 



(76) 



(77) 



(78) 



(79) 



Finally, the comnmtc among themselves. 

As a curiosity, we can actually solve for the operators. The commutation rule suggests = tpj +ig (pi) 
. We now have 



{p2\a\pi) = ^/plSp2+l,p^ = {p2\e '^\pi) 

N 



P2 



-i<p+gip)]/N 



Pi 



J DifiDp exp|i J dt — <^ ^ 



'dp 



I ^ + 1 ) - «5 (P) 



(80) 



where we integrate over paths with p (0) = pi, p (1) = p2- From the integration over ip, we get dp/dt = — 1, 
and so the functional integral vanishes unless p2 = pi — 1, as expected. Finally 



so, if G (p) is the primitive of g (p) , 



so 



, f , dp g{p) 



G{p) = G{p-l) + -\n[p] 



G(p) = -lnr[p+l] 



(81) 



(82) 



(83) 



and 



9{p) = ^i^[p + M 



(84) 



If p ^ 1, we have the Stirling approximation InP [x] ~ (a; — 1) In (x — 1) — (a; — 1) and so ^ [x] ~ In (x — 1) , 
as expected. If we reinstate the h factors, we find we must replace p by p/h, so in the semiclassical limit we 
always have p^ 1. 
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Let us close this Section with a word on the path of integration and measure appropriate to the path 

integral representation. Since the {Ci,Pi) arc canonically conjugated variables, the measure of integration is 
the Liouville measure at each time slice in the path integral, d^i A dpi. To reduce this to a more familiar 
form, we observe that since g (pi) — (— i) {^i — ^*) /2, then d^i A dpi ^ i d^i A d^* /2g' {pi) . Therefore at each 
time slice we must integrate over the whole complex plane; if we adopt the noncanonical (but more usual) 
(^i,^*) pair as independent variables, then a non trivial measure arises. 

This subtlety will not be an obstacle in what follows, since the relevant expectation values will be com- 
puted directly from symmetry arguments or by using the properties of the Heisenberg (g-number) operators 
involved, rather than by an explicit evaluation of the path integral. 

B. The dynamics in the new variables 

At this point we introduce the eigenvectors fpj of the matrix Jij and the corresponding eigenvalues jp. 
For example, consider the case in d = 1, Ng = 2K + 1. 



(85) 



The eigenvectors of Jy are 



exp 



2TTipj 



-K<p<K 



(86) 



and the eigenvalues 



jp = 2 J cos 



2Trp 



(87) 



We now split all variables into a homogeneous and an inhomogeneous part 



(88) 



(89) 



and similarly 



pi 



(90) 
(91) 



The Hamiltonian becomes 



H{pl,rl,X;,) = -plY^Jij [expi [X] - X]]] -Y^Ji^r] [expi [X] - X^]] 



(92) 
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V. THE LINEARIZED APPROXIMATION 



We see from Ea. (|92ll that the model ressembles the dynamics of a sohd, with the Xi playing the role of 
the ion positions and a per iodic interaction potential between ions. There is a sophisticated technology to 
deal with such systems |4l|. The simplest possible approach is the linearized approximation in which the 
excitations of the "solid" are described as a free phonon gas. We shall now develop the implications of this 
view. 

In this Section, we shall derive concrete expressions for the Heisenberg operators corresponding to the Xp 
and Tp operators introduced above, Eqs. (|89|) and H91|) . These expressions shall be used in the next Section 
to derive an analytic approximation for the one body density matrix. 

A. The lowest order equilibrium theory 

We obtain the lowest order theory by keeping only the quadratic terms in the classical action. It is a 
theory of linear fields, and so we may either attempt to solve the equations of motion for the propagators, 
or else solve the Heisenberg equations, which are the same as the classical equations of motion, and compute 
the propagators afterwards. We shall adopt the second path. 

In this Section, therefore, we shall work directly in terms of g-number operators, the Heisenberg equations 
and canonical commutation relations, rather than from the path integral. 

The "free" quadratic part of the Hamiltonian is 

H {po, Tp, Xp) = XI ^ ^^^^'^ ^ ^ ^ "^'^ ~ ~ * ^ "^'^^^ ~ 

i ij ij 

+7V,^Po(po-l) (93) 

Call 



j^p-2(2J-jp) = 8Jsin2^ (94) 



X [X, -X,f^2Y, [(Xjf - X,X.,\ = J2 ^pX-pXp (96) 

ij ij P 

SO (assuming for po the constrained value po = n) 



U 

~2 / y ■ ^p' p ' 2 / ^ ' p — p' 'p 2 
p p p 

The Heisenberg equations of motion are just the classical Hamilton equations 



H {rp, Xp) = ^Y^ r^pTp + ^ X VpX^pXp - VpT-pXp (97) 



(98) 



We seek a solution of the form 



dvp i 



(recall that the Xi are not hermitean, but they satisfy {X*)^ = (X^p)*) 



therefore 



Tp (t) = i-t) 



Cpe-"^"' - Clpe"^"' 



A: 



U 



-C„ 



p — l^p^p 

OJp- — 



with a dispersion relation 



R - ^ r 



Up = y Up (Uji + ^) 
We now write down the canonical commutation relations (cfr. Eq. (j79|l ) 



[n,x,]^{~t) 



whereby 



[rp,Xq] = (-i)(5p+,,o 
Substituting our solution of the Heisenberg equations, we obtain [Cp,Cq] 



Therefore 



nvp 



where is a canonical destruction operator. The final formulae read 
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With the same argument, we get (X*)^ 



L0„ + 



OJr, 



(110) 



COp + ■ 



(111) 



B. Wick theorem and boundary conditions 



Because we have been able to keep the discussion so far at the level of the Hamiltonian and canonical 
operators, we did not have to deal with the issue of the periodicity conditions on the phase variables. 
However, in order to actually use this theory to compute physical expectation values, we must confront this 
issue. 

Let us first ignore all periodicity conditions. This means that we are concerned with the expectation values 
(0) of the covering theory. Because the lowest order approximation to the Hamiltonian is Gaussian, Wick's 
theorem implies 



((e*^)) = e-((^'))/2 (112) 

for any linear function A of the Xp's and Tp's. 

Now we turn to the issue of the expectation values in the physical theory, where the periodicity conditions 
are enforced. To this effect it is convenient to think in terms of the path integral representation of the 
expectation values. 

One way to implement the boundary conditions is not to simply match the history at the second branch 
to the initial state at time i*^^^ = 0, but to all possible translations of this initial state, namely, translating 
all the phases by an amount tpi ^ (pi + 27rfci, for all possible integer ki. Since the translation operator for 
phases is pi, this amounts to defining a new expectation value 

(e*^) =C^((e2"^S.'='^^"(°)e*^)). (113) 

In this equation, the expectation value on the left hand side is physical, and the one on the right hand 

side belongs to the covering theory. Observe that because we are using a CTP representation of the path 
integrals, there are no ordering ambiguities; this is one of the side benefits of the CTP formulation, even 
when discussing equilibrium properties. 

Adding the new exponentials is equivalent to forcing pl (0) to be an integer. Since the homogeneous 

value pIj ^ is already constrained to be an integer by the integration over the Lagrange multiplier p'--^^ (0) , it 
may be ignored. So we shall adopt the definition 

(e^^) = c5^((e^'^^S.'=^'-f'(0)e^^)) (114) 

ki 

The constant C is defined by the condition that (1) = 1, so 



C ^ exp i -277^ J2 kiMimkm \ = 1 



(115) 



ki 



Im 
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where 



Mlm = {{riTm)) (116) 

We now have 



(e"> = Ci:e-(([^«-S.'-i-<"'l)) (U7) 



The exponent reads 



= {{A^))+4nJ2k^ ((rf (0)a)) + (27r)' ^ fc,M;„fc„ (118) 



,.(2) /n^ a\\ _i_ ' 

Im 

Finally, observe that 



M;„ = -J^5^e2-f('— )/^^^ (119) 

p/0 



P 



Since the sum is restricted to nonvanishing momenta, Mim has a zero mode, corresponding to homogeneous 
configurations. In the orthogonal subspace, Mim has an inverse 

j^-l ^ J_ ^2nip{l-m)/N, ^ 

Let us check the asymptotic form of these expressions 

1. The SF limit 

The matrix Mi^ represents the particle number fluctuation correlations between sites I and m. In the SF 
limit, we expect Mi^ = n Si^ — ■ Indeed, in this limit 2up — > Up. So 

(e^^) ~ Ce^ii^')) ^e-'"E^'=^(('^'*''(°)^»+'^"'=''] (121) 

ki 

The sum will be dominated by the fc, = term, so 



In other words, treating the density variables as continuous (thereby ignoring periodicity conditions) is 
not a serious mistake. This is consistent with the large number fluctuations in this regime. 
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2. The MI limit 



For the same reasons, in the MI hmit we expect Mim —>■ 0. Then many terms contribute to the sum, and 
we may replace the discrete sum by an integral. Under this approximation, in fact, we are forcing the not 
just to be an integer, but to vanish. 

Completing the square, we get 




(123) 



As expected, in this limit 




We can prove this by observing that, for A — Pr^, (^yl^^ (0) ~ PMhk. This implies that the particle 
number fluctuations (r^ (0) ri (0)) vanish in this limit. 



VI. THE ONE-BODY DENSITY MATRIX 



We may now turn to computing the one-body density matrix Eq. 



^iik = {aj ih) afe ih)) = (expz [^f* - Cl] (ti)) (125) 

Observe that in our variables, the observable to be computed is a pure exponential: there are no square 
roots to be developed. This is the whole point of introducing the new variables. 

We shall use the machinery introduced above. The desired expectation value is given by Eqs. H117|l and 
(TO 



A. The one body density matrix in the covering theory 



The first step is to compute the expectation value disregarding periodicity conditions, that is, extending 
the path integral to the covering space 



{{4 {h)au (ii))) ^ ((exp^ [Cf - a] [h))) (126) 

To compute this expression, one would aim to split the ^ variables into homogeneous and inhomogeneous 
parts, and to use the formulae above. However, there is a difficulty. We have solved for the Heisenberg 
operators corresponding to the inhomogeneous parts, but not for those of the homogeneous terms. As it 
happens, however, using a symmetry argument saves this effort. Because of the symmetry and the total 
particle number constraint we must have 



o-ioo = {a\ (^i) "fc (^i)) = (127) 

Using this property we can avoid computing explicitly the expectation values for the homogeneous operators. 
Notice however that Eq. (|127() involves a physical expectation value (as oppossed to an expectation value 
within the covering theory). 
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To return to our main argument, we shall seek an expression not for the expectation value 
(ti) Ofc but rather for the ratio between this expectation value and the occupation number 

a|. (ti) ak (^i)^^ • obtain an expression for this ratio, we shall exploit the fact that we are seeking 

the expectation value of an exponential. This expectation value is equivalent to a generating functional for 
connected graphs W where we couple to a source ju (t) — —6ikS [t — ti) and ^f* to j2i (t) = SuSp {t — ti) . 
Formally 

((at(ti)afc(ii)))-e''^[^-^-l (128) 

Decomposing in modes, we see that this is the same as coupling the homogeneous terms to sources 
— =F(5 (i — ti) (independent of k and I) and the inhomogeneous terms to sources ri['^ (t) = —fpkS {t — ti) 
and rii'l it) ^fpiS{t-h). 

In the diagonal case, we would have k = I, and we would find the expectation value on symmetry grounds 
alone. So now we aim to identify how the nondiagonal case is different from the diagonal one. With this goal, 
we write r][p' (t) — r][p (t) + Sri[p'^^ (t) , Sri[p^^ {t) — Sri[p^^6 (t — ti). W has the functional Taylor expansion 



W 



■ .0 (fe) (/) 

J a 1 Vlp ; '?2p 



E ^ E - E Gp....P, ih, -h) svilf.Mif (129) 

g>0 ^' Pl Pq 

where the Gp-^,,,p^ {ti, ...ti) are the time-ordered connected expectation values 

Gp,...p, ih, ...%) = (ti) ...Xl (t,)))^ (130) 

computed for a field driven by the sources Jaj'?ip)^2p- Since the correlation functions are continuous, and 

these sources turn on at the same time as the ^r/^p they may be ignored. Keeping only up to quadratic 
terms, we get 

a\ ih) ak ih))) = e'^« exp | (^^) ^ {{XpX^p)) |/,, - /..j^j (131) 



Wo = W 



' W W 

J a I Vlp I 



(132) 



Observe that we have separated the diagonal and non-diagonal contributions. The former, encoded into 
Wo, shall be obtained below from the symmetry condition Eq. H127f) . so we need not worry about explicitly 
computing the path integral. 



In the vacuum state 



{{XpX,)) = {{X;X*^)) ^ Sp+,^^ (133) 



1/ f |2 4 . 2 TTp (fc - I) 

\.fpk-fpi\ ==]^sm — (134) 
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so 



{la\{h)au{h)))^e''^° X[l-k] (135) 

where 



B. Periodicity corrections 

We now consider the further corrections in Eqs. I|117l) and H118() coming from the periodicity conditions. 
For this, we need to evaluate (^(j'f^ (0) Aki (^i)^^ ^i^h 

Aki=^f*-ek (137) 

First observe that there is no loss of generality in taking ti = 0+. Also there is no problem here with the 
homogeneous terms, because they commute with ap and so give a vanishing expectation value. 
There is no loss of generality in setting the site index A; = 0. In vacuum, 



5oj - S,nj - - [Mj„, + Mjo] 



n 

By symmetry, we must have ctioo = n. We use this condition to determine Wq, getting 



(138) 



Y \m] 

ai„o = nX H (139) 

where X [m] is defined in Eq. H136|l and 

Y[m]^2^e L Z^. ' ' ^ Z^.j (140) 

where vj^ was introduced in Eq. (|138|l . 

These expressions determine the one-body density matrix in the two limiting cases U / J ^ Q and U / J — *■ oo. 
In the former, corresponding to the deep superfluid phase, we have y [m]/y[0] ^ 1 and <^Wm\u/ j^i = 
nX [to] , given by Eq. 1] 136(1 . 

In the latter limit, corresponding to the deep Mott insulator phase, we may replace the sum in Eq. H14U|) 
by an integral over continuous variables ki . In this continuum approximation we get 

crimo|,7/j»i = n (X [m\f (141) 
In the intermediate region, we may approximate Eq. H140|l by 

Y H = n (^C;^,e-2-'^p) ^3 (T,D^,e-^^"^-) (142) 

p>0 
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where i?3 is the Elhptic Theta function [6£ 



2E 



27rp/ 



npm 



Up 2 
— COS 

2ajp 



npm 



(143) 




sm 



2aj„ 



sm 



2Trpm 



(144) 




2Trpl 



(145) 



sm 



2Trpl 



(146) 



This approximation is further discussed in Appendix B. We have checked its accuracy for smah lattices, 
by comparing it to a numerical evaluation of Eq. p4()(l . 



VII. RESULTS AND REMARKS 



A. Results 



In this Section, we shall compare the analytic results above against an exact calculation of the momentum 
distribution function, Eq. 0, for an one dimensional lattice of 9 sites and 9 atoms (n = 1). The exact 
solution was obtained by numerical diagonalization of the Bose Hubbard Haniiltonian. We set J = 1 and 
change U from to 60. We have performed similar calculations for 5 and 7 sites, finding the results to be 
totally consistent with the iV = 9 case. The allowed values of momentum are given by p{q) — q {2TTh/aNs) 
, where a is the lattice spacing and q is an integer. Since by symmetry p{—q) = p{q), there are only five 
independent occupation numbers, corresponding to g = (the condensate) to 4. These are plotted in Figs. 
1 to 5, respectively. 

In these figures we have also plotted the occupation numbers as given by the PNC method ( Bogoliubov) 
calculations , and by first order strong coupling perturbation theory. In all the plots, The solid-red line is 
our prediction, the black dash-dotted line is the exact numerical solution, the pink dots correspond to first 
order strong coupling perturbation theory, and the blue-dashed line to the PNC method. 

We see that for these small lattices our model fares worse than perturbation theory or the PNC approach 
in the corresponding limits of the deep Mott or superfluid regions, but unlike these formalisms, it sustains 
an uniform accuracy throughout. It therefore achieves the goals set in the Introduction. 



B. Remarks 



In this paper we have presented an analytic approximation for the one body density matrix (or equivalently, 
its Fourier transform, the momentum distribution) for a cold gas of structureless bosons in an homogeneous 
optical lattice. We have focused on the regime of low integer filling factor near the superfluid - insulator 
transition, which is not sufficiently covered in the literature. We have checked our results against exact 
calculations for small lattices, and against the theoretical predictions from the Bogoliubov approach and first 
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FIG. 1: The occupation number for the homogeneous mode, as a function of U; n = 1, Ns = 9 and J = 1. The 
soUd-red Une is our prediction, the black dash-dotted line is the exact numerical solution. The pink dots correspond 
to first order strong coupling perturbation theory, and the blue-dashed line to the PNC method. 
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FIG. 2: The occupation number for the first mode as a function of U; the conventions are the same as in Fig. 1. 

order strong coupling perturbation theory. Our model interpolates between these theoretical alternatives, 
keeping an uniform accuracy in the transition region. 

Our model works deep in the MI region, because it is built to be exact when J — 0. This is an advantage 
of our choice of variables over the usual density -phase variables. In the superfluid regime the model predicts 
that quantum fluctuations in the Xp degrees of freedom scale like U (cf. Eq. and thus also qualitatively 

captures the decay of condensate population and the increase of non-condensate atoms. 

However, the agreement is not perfect. The qualitative but not quantitative agreement suggest that 
higher order corrections are required for a proper description of the physics. For observables like the number 
fluctuations at one site, which vanish in the Mott regime according to the linearized approximation, for 
example, higher order corrections would be dominant. 

Quantum corrections will also be important for the dynamic structure factor. There is no contradiction 
between the phonon spectrum of our model (cfr. Eq. (|1U4|| ') and a gapped dynamic structure factor, because 
in the Mott regime the amplitudes of the single phonon poles go to zero, while other poles arise because of 
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FIG. 3: The occupation number for the second mode as a function of U ; the conventions are the same as in Fig. 1. 
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FIG. 4: The occupation number for the third mode as a function of U; the conventions are the same as in Fig. 1. 



higher order corrections. However, in this paper, we have not presented actual results for the particle number 
fluctuations nor the dynamic structure factor; these must be included in the list of unfinished business we 
discussed in the introduction. 

A preliminary comparison we made against available experimental results |60ll6lj of the condensate fraction 
from an array of one-dimensional lattices contained within a three dimensional trap for variable U/ J showed 
fair agreement between the experimental results and the predictions of our model. In these experiments, the 
central tubes had around Ng — 60 populated sites [gJ. The mean occupation number was close to n = 2 
near the center of the trap, and close to n = 1 if averaged over all lattices (T^. We have compared the 
experimental results to the predictions of our model for several values of Ne around 60, and filling fractions 
n = 1 and 2. The results are fairly independent of Ng in this range, and very sensitive to n instead. As a 
typical representative, we show in Fig. 6 the prediction of our model for the condensate fraction for Ng = 61 
and n — l.We have superimposed the experimental results as reported in |60j| . 

We do not regard this as a validation of our model, since it was derived for a translationally invariant 
lattice and the parabolic confinement is not adequately included in our model. Nevertheless, the agreement 
is encouraging and suggest that our model might be more suitable for trapped systems as in this case, in 
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FIG. 5: The occupation number for the fourth mode as a function of U; the conventions are the same as in Fig. 1. 

condensate fraction 




FIG. 6: Condensate fraction (%) plotted against U . Red solid line: prediction from our model using the parameters 
n = 1, A'^a = 61 and J = 1, Dots: Experimental points obtained from Fig. 4a in Ref |60||. 

contrast to the commensurate translationally invariant lattice, there is not a sharp MI transition. We defer 
a detailed discussion to a future communication '62^. 

In summary, in our view, the most important contribution of this work is that it is the first step in 
the formulation of a quantum field theoretical approach capable of dealing with the intermediate regime. 
Even though the agreement with exact numerical solutions is not perfect, we find it satisfactory because 
we are using only the first order approximation. It is a reasonable expectation that by including higher 
order corrections we might narrow the present gap. We are perhaps still a long way from a reliable, fully 
nonequilibrium model of the initialization process of a QIP device based on cold atoms on an optical lattice, 
but from this work we have gained some confidence that we are moving in the right direction. 
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APPENDIX A: APPROXIMATE APPROACHES TO THE MOMENTUM DISTRIBUTION 

FUNCTION 

In this appendix we shall derive the formulae we plotted in the Figures to match against our model. We 
include it only to dispell any ambiguity regarding notation. 

1. Strong coupling Rayleigh - Schrodinger perturbation theory 

This is just ordinary perturbation theory in the parameter J, starting from the state \MI) in Eq. |(7J). 
The BH Hamiltonian Eq. |Q is written as H = Hq + Hi, where Hq is the U term and Hi — — J^ij Jij^\o'j- 
Since MI) = 0, the vacuum energy is unchanged to first order. Hi \MI) is a superposition of one 

particle-hole states, all of which have energy U above the vacuum, so the first order ground state is 



|T> - \MI) - -Hi \MI) 



(Al) 



and the momentum distribution function is 



4^ / 

Uq — n-\- 'jj'TT' (n + I) COS 
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(A2) 



2. The PNC method 

To simplify the problem, we shall consider only the case of a homogeneous, time - independent lattice. 
The starting point of the method is the Heisenberg equation of motion for the destruction operator aj 



[i7,a,(t)]-^J, 



Ua]al 



(A3) 



Parameterize 



-ao 



l + ^e2"*f^'/^=Ap 



(A4) 



There are three key observations: 1) the operators Ap preserve total particle number, 2) the fact that the 
one-body density matrix allows for a homogeneous eigenvector implies ^aJaoAp^ ~ for all p, and 3) we 
have the exact identity 



ajao 



1 + ^AtAp 

p#0 



N 



(A5) 
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Now we develop a perturbative expansion in inverse powers of iV, assuming oq ~ O (^N^ and Ap 
O {l/^/rtj . Multiplying Eq. l|X3|) by aj we get, to first order 



(A6) 



(A7) 



As usual, we seek a solution through a Bogoliubov transformation. Taking into account the commutation 
relations for the Ap we get 



(A8) 



2 _ 2 _ 1 
''p '^p 



(in this simple problem, we may assume the Bogoliubov coefficients are real). We get 



(A9) 



LUp-Un — 



Cp — Unsp — 



(AlO) 



Unc„ 



ujp + Un + — 



Sp = 



from where we recover the dispersion relation Eq. 11041) and 



(All) 



V2 



1 



1/2 



(A12) 



The momentum distribution function is Up = for p 7^ 0, and uq — N — X^p^^o^p ^'-'^ homogeneous 
mode. 



APPENDIX B: APPROXIMATE FORMULA FOR THE ONE-BODY DENSITY MATRIX 

The idea is to evaluate Eq. (|140|l by decomposing the quadratic term in a sum of squares. Of course, one 
possibility is to write 



(Bl) 



The problem is that the requirement that all kj be integer places a highly nontrivial constraint on the kp. 
Consider instead the functions 



/p (j) = sign 



cos 



(B2) 
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gp U) = sign 



sm 



(B3) 



These functions are not orthogonal, but they are a basis. Therefore we can always write 

kj = «o + ^ X! ("P'^P + ^Pfp 



p>0 



Observe that fo is always a null eigenvector of Mij. We expect the / and g's will be approximate eigen- 
vectors. For a large number of sites, 6^^^**^/^= and fp,gp will be nearly orthogonal unless q = ±p, and we 
shall have 



/p(^)-]^{i + 2E 



cos 



2ttpI 



cos 



(B5) 



So 



where of course = and 



9p I 







'2T:pl' 






2TTpj' 


sin 


[ Ns \ 




^ sin 


[ Ns \ 



I 

^ Mjigp {I) - Apgp [j) 



(B6) 



(B7) 



(B8) 



Then from the decomposition Eq. (jB4|l we get 



(B9) 



i p>0 



(BIO) 



ij p > 



(Bll) 



where the coefficients are given in Eqs. (|143ll . (I144II . (|145|l and (|146|l . 

Although each term of the series Eq. H140|l factorizes, there are correlations among the a and b coefficients 
from the discreteness of the kj. For example, for three sites we have that fei = fci — must be an integer, 
but ai = kQ — ki + (bi/2) will be integer if bi is even or half-integer if bi is odd. However, when the number of 
sites is large we may neglect these correlations, and assume that the a and b coefficients simply take integer 
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values. Under this approximation, the multiple sum Eq. (|14U|I factorizes, and we obtain Eq. (|142|1 . where 
i?3 is the Elliptic Theta function [6^ 

t?3 (z, g) = 1 + 2 ^ g"' cos [2nz] (B12) 

n=l 
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